Ladislas Nalborczyk
LPC, LNC, CNRS, Aix-Marseille Univ.
Cours n°01 : Introduction à l'inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux
Cours n°09 : Modèles multi-niveaux généralisés
Cours n°10 : Data Hackathon
Le modèle linéaire Gaussien qu'on a vu aux Cours n°03 et n°04 est caractérisé par un ensemble de postulats, entre autres choses :
\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma)\\ \mu_{i} &= \alpha + \beta_{1} \cdot x_{1i} + \beta_{2} \cdot x_{2i}\\ \end{align} \]
Cette modélisation (le choix d'une distribution Normale) induit plusieurs contraintes, par exemple :
Comment modéliser certaines données qui ne respectent pas ces contraintes ? Par exemple, la proportion de bonnes réponses à un test (bornée entre 0 et 1), un temps de réponse (restreint à des valeurs positives et souvent distribué de manière approximativement log-normale), un nombre d'accidents…
Nous avons déjà rencontré un modèle différent : le modèle Beta-Binomial (cf. Cours n°02).
\[ \begin{align} w &\sim \mathrm{Binomial}(n, p) \\ p &\sim \mathrm{Beta}(a, b) \end{align} \]
Cette modélisation induit deux contraintes :
Comment pourrait-on ajouter des prédicteurs à ce modèle ?
\[ \begin{align} y_{i} &\sim \mathrm{Binomial}(n, p_{i}) \\ f(p_{i}) &= \alpha + \beta \cdot x_{i} \\ \end{align} \]
Objectifs :
Deux changements par rapport au modèle Gaussien :
Les fonctions de lien ont pour tâche de mettre en correspondance l'espace d'un modèle linéaire (non borné) avec l'espace d'un paramètre potentiellement borné comme une probabilité, définie sur l'intervalle \( [0, 1] \).
Les fonctions de lien ont pour tâche de mettre en correspondance l'espace d'un modèle linéaire (non borné) avec l'espace d'un paramètre potentiellement borné comme une probabilité, définie sur l'intervalle \( [0, 1] \).
La fonction logit du GLM binomial (on parle de “log-odds”) :
\[ \text{logit}(p_{i}) = \log \left(\frac{p_{i}}{1 - p_{i}} \right) \]
La cote d'un événement (odds en anglais) est le ratio entre la probabilité que l'événement se produise et la probabilité qu'il ne se produise pas. Le logarithme de cette cote est prédit par un modèle linéaire.
\[ \log \left(\frac{p_{i}}{1 - p_{i}} \right) = \alpha + \beta \cdot x_{i} \]
Pour retrouver la probabilité d'un événement, il faut utiliser la fonction de lien inverse, la fonction logistique (ou logit-inverse) :
\[ p_{i} = \frac{\exp(\alpha + \beta \cdot x_{i})}{1 + \exp(\alpha + \beta \cdot x_{i})} \]
Ces fonctions de lien posent des problèmes d'interprétation : Un changement d'une unité sur un prédicteur n'a plus un effet constant sur la probabilité mais la modifie plus ou moins en fonction de son éloignement à l'origine. Quand \( x = 0 \), une augmentation d'une demi-unité (i.e., \( x = 0.5 \)) se traduit par une augmentation de la probabilité de \( 0.25 \). Puis, chaque augmentation d'une demi-unité se traduit par une augmentation de \( p \) de plus en plus petite…
Deuxième complication : cette fonction de lien force chaque prédicteur à interagir avec lui même et à interagir avec TOUS les autres prédicteurs, même si les interactions ne sont pas explicites…
Dans un modèle Gaussien, le taux de changement de \( y \) en fonction de \( x \) est donné par \( \partial(\alpha + \beta x)~/~\partial x = \beta \) et ne dépend pas de \( x \) (i.e., \( \beta \) est constant).
Dans un GLM binomial (avec une fonction de lien logit), la probabilité d'un événement est donnée par la fonction logistique :
\[ p_{i} = \frac{\exp(\alpha + \beta \cdot x_{i})}{1 + \exp(\alpha + \beta \cdot x_{i})} \]
Et le taux de changement de \( p \) en fonction du prédicteur \( x \) est donné par :
\[ \frac{\partial p}{\partial x} = \frac{\beta}{2(1 + cosh(\alpha + \beta \cdot x)} \]
On voit que la variation sur \( p \) due au prédicteur \( x \) est fonction du prédicteur \( x \)… !
library(tidyverse)
library(rethinking)
data(chimpanzees) # see ?chimpanzees for more information on the dataset
df1 <- chimpanzees
str(df1)
'data.frame': 504 obs. of 8 variables:
$ actor : int 1 1 1 1 1 1 1 1 1 1 ...
$ recipient : int NA NA NA NA NA NA NA NA NA NA ...
$ condition : int 0 0 0 0 0 0 0 0 0 0 ...
$ block : int 1 1 1 1 1 1 2 2 2 2 ...
$ trial : int 2 4 6 8 10 12 14 16 18 20 ...
$ prosoc_left : int 0 0 1 0 1 1 1 1 0 0 ...
$ chose_prosoc: int 1 0 0 1 1 1 0 0 1 1 ...
$ pulled_left : int 0 1 0 0 1 1 0 0 0 0 ...
On cherche à savoir si la présence d'un singe partenaire incite le chimpanzé à appuyer sur le levier prosocial, c'est à dire l'option qui donne de la nourriture aux deux individus. Autrement dit, est-ce qu'il existe une interaction entre l'effet de la latéralité et l'effet de la présence d'un autre chimpanzé sur la probabilité d'actionner le levier gauche.
Observations (pulled_left) : Ce sont des variables de Bernoulli. Elles prennent comme valeur 0/1.
Prédicteur (prosoc_left) : Est-ce que les deux plats sont sur la gauche ou sur la droite ?
Prédicteur (condition) : Est-ce qu'un partenaire est présent ?
\[ \begin{align} L_{i} &\sim \mathrm{Binomial}(1, p_{i}) \\ \text{(equivalent to)} \quad L_{i} &\sim \mathrm{Bernoulli}(p_{i}) \\ \text{logit}(p_{i}) &= \alpha \\ \alpha &\sim \mathrm{Normal}(0, \omega) \\ \end{align} \]
Modèle mathématique sans prédicteur. Comment choisir une valeur pour \( \omega \)… ?
On écrit le modèle précédent avec brms et on échantillonne à partir du prior afin de vérifier que les prédictions du modèle (sur la base du prior seul) correspondent à nos attentes.
library(brms)
mod1.1 <- brm(
formula = pulled_left | trials(1) ~ 1,
family = binomial,
prior = set_prior("normal(0, 10)", class = "Intercept"),
data = df1,
# stores prior samples
sample_prior = "yes"
)
# extracts prior samples
prior_samples(mod1.1) %>%
# applies the inverse link function
mutate(p = brms::inv_logit_scaled(Intercept) ) %>%
ggplot(aes(x = p) ) +
geom_density(fill = "steelblue", adjust = 0.1) +
labs(x = "Probabilité a priori de tirer le levier gauche", y = "Densité de probabilité")
L'intercept s'interprète dans l'espace des log-odds… pour l'interpréter comme une probabilité, il faut appliquer la fonction de lien inverse. On peut utiliser la fonction rethinking::logistic() ou la fonction plogis().
fixed_effects <- fixef(mod1.2) # effets fixes (ou constants)
rethinking::logistic(fixed_effects) # fonction de lien inverse
Estimate Est.Error Q2.5 Q97.5
Intercept 0.5780701 0.5228402 0.5319469 0.6207176
plogis(fixed_effects) # fonction de lien inverse
Estimate Est.Error Q2.5 Q97.5
Intercept 0.5780701 0.5228402 0.5319469 0.6207176
En moyenne (sans considérer les prédicteurs), il semblerait que les singes aient plus tendance à appuyer sur le levier gauche que sur le levier droit…
library(BEST)
post <- posterior_samples(mod1.2)
intercept_samples <- plogis(post$b_Intercept)
plotPost(intercept_samples, compVal = 0.5, xlab = "Probability of pulling left")
Et si on ajoute des prédicteurs… comment choisir une valeur pour \( \omega \) ?
\[ \begin{align} L_{i} &\sim \mathrm{Binomial}(1, p_{i}) \\ \text{logit}(p_{i}) &= \alpha + \beta_{P} P_{i} + \beta_{C} C_{i} + \beta_{PC} P_{i} C_{i} \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \beta_{P}, \beta_{C}, \beta_{PC} &\sim \mathrm{Normal}(0, \omega) \\ \end{align} \]
pulled_left)# recoding predictors
df1 <- df1 %>%
mutate(
prosoc_left = ifelse(prosoc_left == 1, 0.5, -0.5),
condition = ifelse(condition == 1, 0.5, -0.5)
)
priors <- c(
set_prior("normal(0, 1)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b")
)
mod2.1 <- brm(
formula = pulled_left | trials(1) ~ 1 + prosoc_left * condition,
family = binomial,
prior = priors,
data = df1,
sample_prior = "yes"
)
prior_samples(mod2.1) %>%
mutate(
condition1 = plogis(Intercept - 0.5 * b),
condition2 = plogis(Intercept + 0.5 * b)
) %>%
ggplot(aes(x = condition2 - condition1) ) +
geom_density(fill = "steelblue", adjust = 0.1) +
labs(
x = "Différence dans la probabilité a priori de tirer le levier gauche entre conditions",
y = "Densité de probabilité"
)
priors <- c(
set_prior("normal(0, 1)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b")
)
mod2.2 <- brm(
formula = pulled_left | trials(1) ~ 1 + prosoc_left * condition,
family = binomial,
prior = priors,
data = df1,
sample_prior = "yes"
)
prior_samples(mod2.2) %>%
mutate(
condition1 = plogis(Intercept - 0.5 * b),
condition2 = plogis(Intercept + 0.5 * b)
) %>%
ggplot(aes(x = condition2 - condition1) ) +
geom_density(fill = "steelblue", adjust = 0.1) +
labs(
x = "Différence dans la probabilité a priori de tirer le levier gauche entre conditions",
y = "Densité de probabilité"
)
summary(mod2.2)
Family: binomial
Links: mu = logit
Formula: pulled_left | trials(1) ~ 1 + prosoc_left * condition
Data: df1 (Number of observations: 504)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 0.32 0.09 0.15 0.50 1.00 4549
prosoc_left 0.54 0.18 0.21 0.88 1.00 4499
condition -0.19 0.18 -0.55 0.17 1.00 4438
prosoc_left:condition 0.16 0.35 -0.52 0.86 1.00 4665
Tail_ESS
Intercept 3080
prosoc_left 3135
condition 3172
prosoc_left:condition 3123
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Le modèle linéaire ne prédit pas directement la probabilité mais le log-odds de la probabilité :
\[ \begin{align} \text{logit}(p_{i}) &= \log \left(\frac{p_{i}}{1 - p_{i}} \right) = \alpha + \beta x_{i} \\ \end{align} \]
On peut donc distinguer et interpréter deux types d'effets.
Effet relatif : L'effet relatif porte sur le logarithme du rapport des probabilités. Il indique une proportion de changement induit par le prédicteur sur les chances de succès (ou plutôt, sur la cote). Cela ne nous dit rien de la probabilité de l'événement, dans l'absolu.
Effet absolu : Effet qui porte directement sur la probabilité d'un événement. Il dépend de tous les paramètres du modèle et nous donne l'impact effectif d'un changement d'une unité d'un prédicteur (dans l'espace des probabilités).
Il s'agit d'une proportion de changement induit par le prédicteur sur le rapport des chances ou “cote” (odds). Illustration avec un modèle sans interaction.
\[ \begin{align} \log\left(\frac{p_{i}}{1 - p_{i}}\right) &= \alpha + \beta x_{i} \\ \frac{p_{i}}{1 - p_{i}} &= \exp(\alpha + \beta x_{i}) \\ \end{align} \]
La cote proportionnelle \( q \) d'un événement est le nombre par lequel la cote est multipliée lorsque \( x_{i} \) augmente d'une unité.
\[ \begin{align} q = \frac{\exp(\alpha + \beta(x_{i} + 1))}{\exp(\alpha + \beta x_{i})} = \frac{\exp(\alpha) \exp(\beta x_{i}) \exp(\beta)}{\exp(\alpha) \exp(\beta x_{i})} = \exp(\beta) \end{align} \]
Lorsque \( q = 2 \) (par exemple), une augmentation de \( x_{i} \) d'une unité génère un doublement de la cote.
L'effet relatif d'un paramètre dépend également des autres paramètres. Dans le modèle précédent, le prédicteur prosoc_left augmente le log de la cote d'environ 0.54, ce qui se traduit par une augmentation de la cote de \( \exp(0.54) \approx 1.72 \) soit une augmentation d'environ 72% de la cote.
Supposons que l'intercept \( \alpha = 4 \).
prosoc_left, on obtient \( \text{logit}^{-1}(4 + 0.54) \approx 0.99 \).Une augmentation de 72% sur le log-odds se traduit par une augmentation de seulement 1% sur la probabilité effective… Les effets relatifs peuvent conduire à de mauvaises interprétations lorsqu'on ne considère pas l'échelle de la variable mesurée.
fixef(mod2.2)
Estimate Est.Error Q2.5 Q97.5
Intercept 0.3246441 0.09133685 0.1490843 0.5015520
prosoc_left 0.5408226 0.17914704 0.2060616 0.8831652
condition -0.1903554 0.18019174 -0.5501594 0.1679407
prosoc_left:condition 0.1620778 0.34586825 -0.5169865 0.8597904
post <- posterior_samples(mod2.2)
plotPost(exp(post$b_prosoc_left), compVal = 1, xlab = "Odds ratio")
L'effet absolu dépend de tous les paramètres du modèle et nous donne l'impact effectif d'un changement d'une unité d'un prédicteur (dans l'espace des probabilités).
model_predictions <- fitted(mod2.2) %>%
data.frame() %>%
bind_cols(df1) %>%
mutate(condition = factor(condition), prosoc_left = factor(prosoc_left) )
Ces données représentent le nombre de candidatures à l'université de Berkeley par sexe et par département. Chaque candidature est acceptée ou rejetée et les résultats sont agrégés par département et par sexe.
library(rethinking)
data(UCBadmit)
(df2 <- UCBadmit)
dept applicant.gender admit reject applications
1 A male 512 313 825
2 A female 89 19 108
3 B male 353 207 560
4 B female 17 8 25
5 C male 120 205 325
6 C female 202 391 593
7 D male 138 279 417
8 D female 131 244 375
9 E male 53 138 191
10 E female 94 299 393
11 F male 22 351 373
12 F female 24 317 341
Existe-t-il un biais de recrutement lié au sexe ?
On va construire un modèle de la décision d'admission en prenant comme prédicteur le sexe du candidat.
\[ \begin{align} \text{admit}_{i} &\sim \mathrm{Binomial}(n_{i}, p_{i}) \\ \text{logit}(p_i) &= \alpha + \beta_{m} m_{i} \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \beta_{m} &\sim \mathrm{Normal}(0, 1) \\ \end{align} \]
Les variables :
admit)applications)1 = male)priors <- c(set_prior("normal(0, 1)", class = "Intercept") )
mod3 <- brm(
formula = admit | trials(applications) ~ 1,
family = binomial(link = "logit"),
prior = priors,
data = df2,
sample_prior = "yes"
)
priors <- c(
set_prior("normal(0, 1)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b")
)
# dummy-coding
df2$male <- ifelse(df2$applicant.gender == "male", 1, 0)
mod4 <- brm(
formula = admit | trials(applications) ~ 1 + male,
family = binomial(link = "logit"),
prior = priors,
data = df2,
sample_prior = "yes"
)
summary(mod4)
Family: binomial
Links: mu = logit
Formula: admit | trials(applications) ~ 1 + male
Data: df2 (Number of observations: 12)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.83 0.05 -0.93 -0.73 1.00 2180 2653
male 0.61 0.06 0.48 0.73 1.00 2622 2740
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Être un homme semble être un avantage… ! Le rapport des cotes est de \( \exp(0.61) \approx 1.84 \).
Calculons la différence de probabilité d'admission entre hommes et femmes.
post <- posterior_samples(mod4)
p.admit.male <- plogis(post$b_Intercept + post$b_male)
p.admit.female <- plogis(post$b_Intercept)
diff.admit <- p.admit.male - p.admit.female
plotPost(diff.admit, compVal = 0)
On examine les prédictions du modèle par département.
Les prédictions du modèle sont très mauvaises… Il n'y a que deux départements pour lesquels les femmes ont de moins bonnes prédictions que les hommes (C et E) alors que le modèle prédit une probabilité d'admission plus basse pour tous les départements…
Le problème est double :
C'est le paradoxe de Simpson… remarques :
On construit donc un modèle de la décision d'admission en fonction du genre, au sein de chaque département.
\[ \begin{align} \text{admit}_{i} &\sim \mathrm{Binomial}(n_{i}, p_{i}) \\ \text{logit}(p_i) &= \alpha_{\text{dept}[i]} + \beta_{m} m_{i} \\ \alpha_{\text{dept}[i]} &\sim \mathrm{Normal}(0, 1) \\ \beta_{m} &\sim \mathrm{Normal}(0, 1) \\ \end{align} \]
# modèle sans prédicteur
mod5 <- brm(
admit | trials(applications) ~ 0 + dept,
family = binomial(link = "logit"),
prior = set_prior("normal(0, 1)", class = "b"),
data = df2
)
# modèle avec prédicteur
mod6 <- brm(
admit | trials(applications) ~ 0 + dept + male,
family = binomial(link = "logit"),
prior = set_prior("normal(0, 1)", class = "b"),
data = df2
)
summary(mod6)
Family: binomial
Links: mu = logit
Formula: admit | trials(applications) ~ 0 + dept + male
Data: df2 (Number of observations: 12)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
deptA 0.68 0.10 0.50 0.88 1.00 2362 2872
deptB 0.64 0.11 0.42 0.86 1.00 2460 2719
deptC -0.58 0.07 -0.73 -0.43 1.00 3659 2532
deptD -0.61 0.09 -0.79 -0.44 1.00 2941 2615
deptE -1.05 0.10 -1.24 -0.86 1.00 4159 2991
deptF -2.57 0.15 -2.87 -2.28 1.00 4039 2815
male -0.10 0.08 -0.26 0.05 1.00 1894 2583
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
fixef(mod6)
Estimate Est.Error Q2.5 Q97.5
deptA 0.6839440 0.09729755 0.4970922 0.87843033
deptB 0.6384439 0.11247404 0.4165107 0.85970950
deptC -0.5775175 0.07464361 -0.7268685 -0.42740448
deptD -0.6097237 0.08621418 -0.7880469 -0.44189251
deptE -1.0513581 0.09739745 -1.2433975 -0.86129138
deptF -2.5722057 0.14858430 -2.8693439 -2.28378476
male -0.1034019 0.07873304 -0.2563696 0.05302116
Maintenant, la prédiction pour \( \beta_{m} \) va dans l'autre sens… La rapport des cotes (odds ratio) est de \( \exp(-0.1) \) soit 90% de la cote des femmes.
Les hommes et les femmes ne postulent pas aux mêmes départements et les départements varient par leur probabilité d'admission. En l'occurrence, les femmes ont plus postulé aux départements E et F (avec une probabilité d'admission plus faible) et ont moins postulé aux départements A ou B, avec une probabilité d'admission plus haute.
Pour évaluer l'effet du sexe sur la probabilité d'admission, il faut donc se poser la question suivante : “Quelle est la différence de probabilité d'admission entre hommes et femmes au sein de chaque département ?” (plutôt que de manière générale).
Retenir que le modèle de régression peut être généralisé à différents modèles de génération des données (i.e., différentes distributions de probabilité, comme la distribution Normale, Binomiale, Poisson, etc) et que l'espace des paramètres peut être “connecté” à l'espace des prédicteurs (variables mesurées) grâce à des fonctions de lien (e.g., la fonction logarithme, exponentielle, logit, etc).
Retenir la distinction entre effet relatif (e.g., un changement de cote) et effet absolu (e.g., une différence de probabilité).
Travailler avec des sujets humains implique un minimum de coopération réciproque. Mais ce n'est pas toujours le cas. Une partie non-négligeable des étudiants qui s'inscrivent pour passer des expériences de Psychologie ne se présentent pas le jour prévu… On a voulu estimer la probabilité de présence d'un étudiant inscrit en fonction de l'envoi (ou non) d'un mail de rappel (cet exemple est présenté en détails dans deux blogposts, accessibles ici, et ici).
df3 <- read.csv("data/absence.csv")
df3 %>% sample_frac %>% head(10)
day inscription reminder absence presence total
1 Wednesday doodle no 6 11 17
2 Tuesday doodle yes 1 7 8
3 Friday doodle no 7 11 18
4 Thursday doodle no 3 11 14
5 Wednesday doodle yes 0 4 4
6 Friday doodle yes 0 2 2
7 Tuesday doodle no 4 10 14
8 Friday panel yes 0 10 10
9 Monday doodle no 5 4 9
10 Monday panel yes 6 12 18
Écrire le modèle qui prédit la présence d'un participant sans prédicteur.
\[ \begin{aligned} y_{i} &\sim \mathrm{Binomial}(n_{i}, p_{i}) \\ \text{logit}(p_{i}) &= \alpha \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \end{aligned} \]
mod7 <- brm(
presence | trials(total) ~ 1,
family = binomial(link = "logit"),
prior = set_prior("normal(0, 1)", class = "Intercept"),
data = df3,
cores = parallel::detectCores()
)
fixef(mod7) # effet relatif (log de la cote)
Estimate Est.Error Q2.5 Q97.5
Intercept 1.162342 0.1923876 0.7876458 1.542933
fixef(mod7) %>% plogis # effet absolu (probabilité de présence)
Estimate Est.Error Q2.5 Q97.5
Intercept 0.761758 0.5479491 0.6873256 0.8238907
On commence par re-coder en dummy variables le reminder et l'inscription.
df3 <-
df3 %>%
mutate(
reminder = ifelse(reminder == "no", 0, 1),
inscription = ifelse(inscription == "panel", 0, 1)
)
head(df3, n = 10)
day inscription reminder absence presence total
1 Friday 1 0 7 11 18
2 Friday 1 1 0 2 2
3 Friday 0 1 0 10 10
4 Monday 1 0 5 4 9
5 Monday 1 1 2 6 8
6 Monday 0 1 6 12 18
7 Thursday 1 0 3 11 14
8 Tuesday 1 0 4 10 14
9 Tuesday 1 1 1 7 8
10 Tuesday 0 1 0 9 9
Écrire le modèle qui prédit la présence en fonction du rappel.
\[ \begin{aligned} y_{i} &\sim \mathrm{Binomial}(n_{i}, p_{i}) \\ \text{logit}(p_{i}) &= \alpha + \beta \times \text{reminder}_{i} \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \beta &\sim \mathrm{Normal}(0, 1) \\ \end{aligned} \]
Écrire le modèle qui prédit la présence en fonction du rappel.
priors <- c(
set_prior("normal(0, 1)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b")
)
mod8 <- brm(
presence | trials(total) ~ 1 + reminder,
family = binomial(link = "logit"),
prior = priors,
data = df3,
cores = parallel::detectCores()
)
Quel est l'effet relatif du mail de rappel ?
exp(fixef(mod8)[2]) # odds ratio between no-reminder and reminder
[1] 3.004387
Envoyer un rappel augmente proportionnellement les chances de présence (i.e., augmente la cote) par environ \( 3 \).
Quel est l'effet absolu du mail de rappel ?
post <- posterior_samples(mod8) # extracting posterior samples
p.no <- plogis(post$b_Intercept) # mean probability of presence when no reminder
p.yes <- plogis(post$b_Intercept + post$b_reminder) # mean probability of presence when reminder
plotPost(p.yes - p.no, compVal = 0, showMode = TRUE) # plotting it
library(tidybayes)
library(modelr)
df3 %>%
group_by(total) %>%
data_grid(reminder = seq_range(reminder, n = 1e2) ) %>%
add_fitted_draws(mod8, newdata = ., n = 100, scale = "linear") %>%
mutate(estimate = plogis(.value) ) %>%
group_by(reminder, .draw) %>%
summarise(estimate = mean(estimate) ) %>%
ggplot(aes(x = reminder, y = estimate, group = .draw) ) +
geom_hline(yintercept = 0.5, lty = 2) +
geom_line(aes(y = estimate, group = .draw), size = 0.5, alpha = 0.1) +
ylim(0, 1) +
labs(x = "Mail de rappel", y = "Probabilité de présence")
Écrire le modèle qui prédit la présence en fonction du mode d'inscription.
\[ \begin{aligned} y_{i} &\sim \mathrm{Binomial}(n_{i}, p_{i}) \\ \text{logit}(p_{i}) &= \alpha + \beta \times \text{inscription}_{i} \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \beta &\sim \mathrm{Normal}(0, 1) \\ \end{aligned} \]
priors <- c(
set_prior("normal(0, 1)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b")
)
mod9 <- brm(
presence | trials(total) ~ 1 + inscription,
family = binomial(link = "logit"),
prior = priors,
data = df3,
cores = parallel::detectCores()
)
post <- posterior_samples(mod9)
p.panel <- plogis(post$b_Intercept) # mean probability of presence for panel
p.doodle <- plogis(post$b_Intercept + post$b_inscription) # mean probability of presence for doodle
plotPost(p.panel - p.doodle, compVal = 0, showMode = TRUE) # plotting it
La probabilité de présence est augmentée de \( 0.17 \) lorsque l'on s'inscrit sur un panel comparativement à une inscription sur un doodle (effet légèrement plus faible que pour le rappel).
Écrire le modèle complet.
\[ \begin{aligned} y_{i} &\sim \mathrm{Binomial}(n_{i}, p_{i}) \\ \text{logit}(p_{i}) &= \alpha + \beta_{1} \times \text{reminder}_{i} + \beta_{2} \times \text{inscription}_{i} \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \beta_{1}, \beta_{2} &\sim \mathrm{Normal}(0, 1) \\ \end{aligned} \]
priors <- c(
set_prior("normal(0, 1)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b")
)
mod10 <- brm(
presence | trials(total) ~ 1 + reminder + inscription,
family = binomial(link = "logit"),
prior = priors,
data = df3,
cores = parallel::detectCores()
)
summary(mod10)
Family: binomial
Links: mu = logit
Formula: presence | trials(total) ~ 1 + reminder + inscription
Data: df3 (Number of observations: 13)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.00 0.59 -0.11 2.16 1.00 2000 2258
reminder 0.92 0.50 -0.08 1.87 1.00 2250 2437
inscription -0.34 0.56 -1.43 0.72 1.00 2130 2425
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Le mode d'inscription et le mail de rappel semblent avoir moins d'effet dans le modèle complet que dans les modèles simples… pourquoi ?
fixef(mod8) %>% exp
Estimate Est.Error Q2.5 Q97.5
Intercept 1.959509 1.272208 1.223267 3.121618
reminder 3.004387 1.468252 1.401986 6.399431
fixef(mod9) %>% exp
Estimate Est.Error Q2.5 Q97.5
Intercept 6.2546164 1.463510 3.0853122 13.93595
inscription 0.3839219 1.540663 0.1577284 0.88431
fixef(mod10) %>% exp
Estimate Est.Error Q2.5 Q97.5
Intercept 2.7227728 1.805266 0.8949859 8.700254
reminder 2.4974062 1.643716 0.9217869 6.470679
inscription 0.7144684 1.743341 0.2385755 2.045819
On a déjà rencontré ce cas de figure (cf. Cours n°04). Lorsque deux prédicteurs contiennent une part d'information commune, l'estimation des pentes est corrêlée…
posterior_samples(mod10) %>%
ggplot(aes(b_reminder, b_inscription) ) +
geom_point(size = 3, pch = 21, alpha = 0.8, color = "white", fill = "black")
En effet, les données ont été collectées par deux expérimentateurs. L'un d'entre eux a recruté tous ses participants via doodle, et n'envoyait pas souvent de mail de rappel. Le deuxième expérimentateur a recruté tous ses participants via un panneau physique présent dans le laboratoire et envoyait systématiquement un mail de rappel. Autrement dit, ces deux variables sont presque parfaitement confondues.
read.csv("data/absence.csv") %>%
sample_frac %>%
group_by(inscription, reminder) %>%
summarise(n = sum(total) ) %>%
spread(key = reminder, value = n) %>%
data.frame
inscription no yes
1 doodle 72 22
2 panel NA 51